Считываем данные, строим исходный ряд и его периодограмму.
ts_data <- read.csv("MRTSSM4453USN.csv", header = TRUE, as.is = FALSE, sep = ',')[,1:2]
ts_data[,1] <- c(0:359)
ts1 <- ts_data[1:336,]
ts <- ts(ts_data[,2], start = c(1992, 1), frequency = 12)
ts1 <- ts(ts1[,2], start = c(1992, 1), frequency = 12)
ts2 <- ts(ts1[c(1:54)], start = c(1992, 1), frequency = 12)
plot.ts(ts1, type = 'l')
spectrum(ts1, log='no', method='pgram', detrend=FALSE)
spectrum(ts1, log='no', method='pgram')
На периодограмме все частоты имеют значительный вклад, в том числе присутствует частота 1/2 (6/12), то есть, возможно, мы встретим “пилу”.
Также продемонструруем периодограмму белого шума.
data <- rnorm(1000, 0, 1)
plot.ts(data, type="l")
spectrum(data, log='no')
spectrum(data, log='no', taper=0)
Построим ряд, его периодограмму и применим фильтр для подавления частот вне диапазона 0.2–0.45.
plot.ts(ts1, type = 'l')
spectrum(ts1, log='no', method='pgram')
ts_data[1:336,]
## DATE MRTSSM4453USN
## 1 0 1509
## 2 1 1541
## 3 2 1597
## 4 3 1675
## 5 4 1822
## 6 5 1775
## 7 6 1912
## 8 7 1862
## 9 8 1770
## 10 9 1882
## 11 10 1831
## 12 11 2511
## 13 12 1614
## 14 13 1529
## 15 14 1678
## 16 15 1713
## 17 16 1796
## 18 17 1792
## 19 18 1950
## 20 19 1777
## 21 20 1707
## 22 21 1757
## 23 22 1782
## 24 23 2443
## 25 24 1548
## 26 25 1505
## 27 26 1714
## 28 27 1757
## 29 28 1830
## 30 29 1857
## 31 30 1981
## 32 31 1858
## 33 32 1823
## 34 33 1806
## 35 34 1845
## 36 35 2577
## 37 36 1555
## 38 37 1501
## 39 38 1725
## 40 39 1699
## 41 40 1807
## 42 41 1863
## 43 42 1886
## 44 43 1861
## 45 44 1845
## 46 45 1788
## 47 46 1879
## 48 47 2598
## 49 48 1679
## 50 49 1652
## 51 50 1837
## 52 51 1798
## 53 52 1957
## 54 53 1958
## 55 54 2034
## 56 55 2062
## 57 56 1781
## 58 57 1860
## 59 58 1992
## 60 59 2547
## 61 60 1706
## 62 61 1621
## 63 62 1853
## 64 63 1817
## 65 64 2060
## 66 65 2002
## 67 66 2098
## 68 67 2079
## 69 68 1892
## 70 69 2050
## 71 70 2082
## 72 71 2821
## 73 72 1846
## 74 73 1768
## 75 74 1894
## 76 75 1963
## 77 76 2140
## 78 77 2059
## 79 78 2209
## 80 79 2118
## 81 80 2031
## 82 81 2163
## 83 82 2154
## 84 83 3037
## 85 84 1866
## 86 85 1808
## 87 86 1986
## 88 87 2099
## 89 88 2210
## 90 89 2145
## 91 90 2339
## 92 91 2140
## 93 92 2126
## 94 93 2219
## 95 94 2273
## 96 95 3265
## 97 96 1920
## 98 97 1976
## 99 98 2190
## 100 99 2132
## 101 100 2357
## 102 101 2413
## 103 102 2463
## 104 103 2422
## 105 104 2358
## 106 105 2352
## 107 106 2549
## 108 107 3375
## 109 108 2109
## 110 109 2052
## 111 110 2327
## 112 111 2231
## 113 112 2470
## 114 113 2526
## 115 114 2483
## 116 115 2518
## 117 116 2316
## 118 117 2409
## 119 118 2638
## 120 119 3542
## 121 120 2114
## 122 121 2109
## 123 122 2366
## 124 123 2300
## 125 124 2569
## 126 125 2486
## 127 126 2568
## 128 127 2595
## 129 128 2297
## 130 129 2401
## 131 130 2601
## 132 131 3488
## 133 132 2121
## 134 133 2046
## 135 134 2273
## 136 135 2333
## 137 136 2576
## 138 137 2433
## 139 138 2611
## 140 139 2660
## 141 140 2461
## 142 141 2641
## 143 142 2660
## 144 143 3654
## 145 144 2293
## 146 145 2219
## 147 146 2398
## 148 147 2553
## 149 148 2685
## 150 149 2643
## 151 150 2867
## 152 151 2622
## 153 152 2618
## 154 153 2727
## 155 154 2763
## 156 155 3801
## 157 156 2219
## 158 157 2316
## 159 158 2530
## 160 159 2640
## 161 160 2709
## 162 161 2783
## 163 162 2924
## 164 163 2791
## 165 164 2784
## 166 165 2801
## 167 166 2933
## 168 167 4137
## 169 168 2424
## 170 169 2519
## 171 170 2753
## 172 171 2791
## 173 172 3017
## 174 173 3055
## 175 174 3117
## 176 175 3024
## 177 176 2997
## 178 177 2913
## 179 178 3137
## 180 179 4269
## 181 180 2569
## 182 181 2603
## 183 182 3005
## 184 183 2867
## 185 184 3262
## 186 185 3364
## 187 186 3322
## 188 187 3292
## 189 188 3057
## 190 189 3087
## 191 190 3297
## 192 191 4403
## 193 192 2675
## 194 193 2806
## 195 194 2989
## 196 195 2997
## 197 196 3420
## 198 197 3279
## 199 198 3517
## 200 199 3472
## 201 200 3151
## 202 201 3351
## 203 202 3386
## 204 203 4461
## 205 204 2913
## 206 205 2781
## 207 206 3024
## 208 207 3130
## 209 208 3467
## 210 209 3307
## 211 210 3555
## 212 211 3399
## 213 212 3263
## 214 213 3425
## 215 214 3356
## 216 215 4625
## 217 216 2878
## 218 217 2916
## 219 218 3214
## 220 219 3310
## 221 220 3467
## 222 221 3438
## 223 222 3657
## 224 223 3454
## 225 224 3365
## 226 225 3497
## 227 226 3524
## 228 227 4681
## 229 228 2888
## 230 229 2984
## 231 230 3249
## 232 231 3363
## 233 232 3471
## 234 233 3551
## 235 234 3740
## 236 235 3576
## 237 236 3517
## 238 237 3515
## 239 238 3646
## 240 239 4892
## 241 240 2995
## 242 241 3202
## 243 242 3550
## 244 243 3409
## 245 244 3786
## 246 245 3816
## 247 246 3733
## 248 247 3752
## 249 248 3503
## 250 249 3626
## 251 250 3869
## 252 251 5124
## 253 252 3143
## 254 253 3212
## 255 254 3603
## 256 255 3464
## 257 256 3916
## 258 257 3776
## 259 258 3994
## 260 259 4056
## 261 260 3588
## 262 261 3741
## 263 262 4007
## 264 263 5147
## 265 264 3333
## 266 265 3261
## 267 266 3596
## 268 267 3643
## 269 268 4096
## 270 269 3966
## 271 270 4166
## 272 271 4139
## 273 272 3736
## 274 273 4003
## 275 274 4012
## 276 275 5444
## 277 276 3486
## 278 277 3397
## 279 278 3761
## 280 279 3768
## 281 280 4222
## 282 281 4104
## 283 282 4409
## 284 283 4140
## 285 284 3955
## 286 285 4145
## 287 286 4135
## 288 287 5634
## 289 288 3488
## 290 289 3642
## 291 290 3907
## 292 291 3966
## 293 292 4242
## 294 293 4307
## 295 294 4572
## 296 295 4307
## 297 296 4260
## 298 297 4261
## 299 298 4488
## 300 299 5812
## 301 300 3578
## 302 301 3606
## 303 302 4074
## 304 303 4077
## 305 304 4456
## 306 305 4482
## 307 306 4598
## 308 307 4452
## 309 308 4346
## 310 309 4343
## 311 310 4638
## 312 311 5972
## 313 312 3792
## 314 313 3792
## 315 314 4436
## 316 315 4143
## 317 316 4702
## 318 317 4740
## 319 318 4761
## 320 319 4697
## 321 320 4416
## 322 321 4555
## 323 322 4926
## 324 323 6128
## 325 324 3933
## 326 325 3916
## 327 326 4445
## 328 327 4358
## 329 328 4861
## 330 329 4769
## 331 330 4993
## 332 331 5017
## 333 332 4454
## 334 333 4676
## 335 334 5057
## 336 335 6326
bandpass_ts <- bandpass(ts_data[1:336,],fhigh=0.45,flow=0.05,win=2,p=.05,detrend = FALSE)
##
## ----- BANDPASS FILTERING STRATIGRAPHIC SERIES-----
## * Number of data points= 336
## * Sample interval= 1
## * Mean value removed= 3008.012
## * Center of bandpass filter = 0.25
## * 269 pos/neg frequency pairs will be bandpassed
ts_b <- ts(bandpass_ts[,2], start = c(1992, 1), frequency = 12)
plot.ts(bandpass_ts[,2], type = 'l')
spectrum(ts_b, log='no', method='pgram')
Полученный ряд не имеет линейного тренда, так как тренд соответствует низким частотам, а также отсутствует гармоника с частотой 1/2 (\((-1)^n\)).
Построим АЧХ для нескольких фильтров:
afc <- function(filter, omega) {
k <- seq_along(filter) - 1
h <- function(o) sum(rev(filter) * exp(-k*1i * o))
abs(sapply(omega, h))
}
freq <- seq(0, pi, 0.001)
filt <- rep(1, 36)
omega <- freq/2/pi
plot(afc(filt, freq) ~ omega, type = 'l')
filt <- c(1,-1)
plot(afc(filt, freq) ~ omega, type = 'l')
Первый для скользящего среднего с окном 36, второй– переход к разностям.
Смоделируем сильно зашумленный ряд и попробуем выделить необходимую частоту.
Применим скользящее среднее к периодограмме шума, для оценки спектральной плотности:
n <- 200
N <- 1000
wnoise <- rnorm(N, 0, 1)
plot.ts(wnoise, type="l")
B <- spectrum(wnoise, log='no', taper=0, plot=FALSE)
plot(B$freq, B$spec, type='l')
lines(B$freq, movavg(B$spec, n, type='t'), type = 'l', col = 'red')
#spectrum(wnoise, kernel = kernel("daniell", c(11,7,3)), log = "no")
#bandpass_ts2 <- bandpass(cbind(c(0:N-1), wnoise),flow=0.2,fhigh=0.5,win=2,p=.5,detrend = TRUE)
#ts_b2 <- ts(bandpass_ts2[,2], start = c(1992, 1), frequency = 12)
#plot.ts(bandpass_ts2[,2], type = 'l')
#spectrum(ts_b2, log='no', method='pgram')
Для белого шума получается константа.
w0 <- wnoise[1]
wnoise <- wnoise[2:N]
cor <- 0.8
rnoise = Reduce(function(prev_v, next_v) cor * prev_v + next_v * sqrt(1 - cor^2), wnoise, w0, accumulate = T)
plot.ts(rnoise, type="l")
B <- spectrum(rnoise, log='no', taper=0, plot=FALSE)
plot(B$freq, B$spec, type='l')
lines(B$freq, movavg(B$spec, n, type='t'), type = 'l', col = 'red')
#bandpass_ts3 <- lowpass(cbind(c(0:N-1), rnoise),fcut=.05,win=2,p=.2,detrend = TRUE)
#ts_b3 <- ts(bandpass_ts3[,2], start = c(1992, 1), frequency = 12)
#plot.ts(bandpass_ts3[,2], type = 'l')
#spectrum(ts_b3, log='no', method='pgram')
Для красного шума убывающая кривая.
Выделение тренда с помощью скользящего среднего с разной длиной окна (берем длину окна, кратную периоду ряда–12):
ma <- function(x, n){filter(x, rep(1 / n, n), sides = 2)}
plot.ts(ts1, type = 'l')
lines(x = time(ts1), y = ma(ts1, 24), type = 'l', col = 'blue')
lines(x = time(ts1), y = ma(ts1, 48), type = 'l', col = 'red')
Недостатком является необходимость достраивать края.
n <- 48
plot.ts(ts1, type = 'l')
lines(x = time(ts1), y = movavg(ts1, n, type='t'), type = 'l', col = 'red')
lines(x = time(ts1), y = ma(ts1, n), type = 'l', col = 'blue')
Применим HP-фильтр с параметром лямбда = 12960.
plot.ts(ts1, type = 'l')
hp_filt <- hpfilter(ts1, freq=12960, type="lambda")$trend
lines(hp_filt, type = 'l', col = 'red')
Результат похож на скользящее среднее, но тренд выделен и на краях. Недостаток- сложность подбора параметра.
time_ts1 <- time(ts1)
model <- nls(ts1 ~ b1*time_ts1+b2, start = list(b1 = 1, b2 = 3))
plot.ts(ts1, type = 'l')
lines(x = time_ts1, y = predict(model, ts1), type = 'l', col = 'red')
Приближение многочленом первой степени плохо описывает поведение на краях.
model <- nls(ts1 ~ b1*time_ts1^2+b2*time_ts1+b3, start = list(b1 = 1, b2 = 3, b3 = 0))
plot.ts(ts1, type = 'l')
lines(x = time_ts1, y = predict(model, ts1), type = 'l', col = 'red')
Приближение многочленом второй степени неплохо описывает тренд.
expR <- lm(log(ts1) ~ time(ts1))
plot.ts(ts1, type = 'l')
lines(x = time_ts1, y = exp(predict(expR)), type = 'l', col = 'red')
Приближение экспонентой неплохо описывает тренд.
Для выделения тренда рассмотрим метод ssa с длиной окна 12.
ssa.result <- ssa(ts1, L=12)
plot(ssa.result, type = 'vectors', idx=1:12)
plot(ssa.result, type = 'paired')
plot(ssa.result, type = 'series')
plot(ssa.result, type = 'wcor')
parestimate(ssa.result)
## $F1
## period rate | Mod Arg | Re Im
## Inf 0.003176 | 1.00318 0.00 | 1.00318 0.00000
##
## $F2
## period rate | Mod Arg | Re Im
## Inf -0.568060 | 0.56662 0.00 | 0.56662 0.00000
##
## $F3
## period rate | Mod Arg | Re Im
## Inf -0.757734 | 0.46873 0.00 | 0.46873 0.00000
##
## $F4
## period rate | Mod Arg | Re Im
## Inf -2.784229 | 0.06178 0.00 | 0.06178 0.00000
##
## $F5
## period rate | Mod Arg | Re Im
## 2.000 -3.250759 | 0.03874 3.14 | -0.03874 0.00000
##
## $F6
## period rate | Mod Arg | Re Im
## 2.000 -0.150974 | 0.85987 3.14 | -0.85987 0.00000
##
## $F7
## period rate | Mod Arg | Re Im
## 2.000 -0.157958 | 0.85389 3.14 | -0.85389 0.00000
##
## $F8
## period rate | Mod Arg | Re Im
## 2.000 -0.555808 | 0.57361 3.14 | -0.57361 0.00000
##
## $F9
## period rate | Mod Arg | Re Im
## 2.000 -0.877564 | 0.41579 3.14 | -0.41579 0.00000
##
## $F10
## period rate | Mod Arg | Re Im
## Inf -0.113651 | 0.89257 0.00 | 0.89257 0.00000
##
## $F11
## period rate | Mod Arg | Re Im
## Inf -0.199074 | 0.81949 0.00 | 0.81949 0.00000
##
## $F12
## period rate | Mod Arg | Re Im
## 2.000 -0.007688 | 0.99234 3.14 | -0.99234 0.00000
ssa.rec.trend <- reconstruct(ssa.result, groups = list(Trend = c(1)))
plot(ssa.rec.trend)
Трендом будем считать только первую компоненту.
plot.ts(ts1, type = 'l')
lines(ssa.rec.trend$Trend, col='red')
Из рассмотренных методов ssa наилучшим образом справляется с выделением тренда.
Добавим еще lowess
plot.ts(ts1, type = 'l')
#lines(lowess(ts1 ~ time(ts1), f=.05, iter=3L), col='green')
lines(lowess(ts1 ~ time(ts1), f=.1, iter=3), col='blue')
#lines(lowess(ts1 ~ time(ts1), f=.15, iter=3L), col=4)
lines(lowess(ts1 ~ time(ts1), f=.2, iter=3), col='red')
Результаты похожи на ssa, сильных отличий при использовании в lowess 10% и 20% точек нет.
dec_ts1 <- stats::decompose(ts1, type = 'multiplicative')
plot(ts1, type = 'l')
lines(dec_ts1$trend, col='red')
plot(dec_ts1$seasonal)
plot(dec_ts1$random)
spectrum(na.omit(dec_ts1$random), log='no', method='pgram', detrend=TRUE)
Полученный тренд неплохо описывает ряд, на периодограмме видно, что в шуме остались гармоники.
Возьмем t.window=13, при общем количестве периодов 28:
stl_ts1 <- stl(log(ts1), s.window = 5, s.degree = 1, t.window = 13, t.jump = 1, t.degree=1)
spectrum(stl_ts1$time.series[,3], log='no', method='pgram', detrend=TRUE)
#acf(stl_ts1$time.series[,3])
#acf(stl_ts1$time.series[,2])
#acf(stl_ts1$time.series[,1])
plot(stl_ts1)
plot(ts1, type = 'l')
lines(exp(stl_ts1$time.series[,2]), col='red')
Полученный тренд неплохо описывает ряд, на периодограмме видно, что в шуме остались гармоники.
Рассмотрим ssa с максимально возможной длинной окна- половина ряда.
ssa.result <- ssa(ts1, L=336/2)
plot(ssa.result, type = 'vectors', idx = 1:20)
plot(ssa.result, type = 'paired', idx = 2:20)
plot(ssa.result, type = 'series', groups = 1:20)
plot(ssa.result, type = 'series', groups = list(1:20))
plot(ssa.result, type = 'wcor', groups = 1:20)
parestimate(ssa.result, group=list(2:12))
## period rate | Mod Arg | Re Im
## 2.999 0.002750 | 1.00275 2.10 | -0.50225 0.86790
## -2.999 0.002750 | 1.00275 -2.10 | -0.50225 -0.86790
## 5.999 0.002663 | 1.00267 1.05 | 0.50123 0.86839
## -5.999 0.002663 | 1.00267 -1.05 | 0.50123 -0.86839
## 2.400 0.002556 | 1.00256 2.62 | -0.86807 0.50157
## -2.400 0.002556 | 1.00256 -2.62 | -0.86807 -0.50157
## 3.999 0.002536 | 1.00254 1.57 | -0.00055 1.00254
## -3.999 0.002536 | 1.00254 -1.57 | -0.00055 -1.00254
## 11.968 0.002470 | 1.00247 0.52 | 0.86747 0.50245
## -11.968 0.002470 | 1.00247 -0.52 | 0.86747 -0.50245
## -2.000 0.002032 | 1.00203 -3.14 | -1.00203 -0.00000
ssa.rec.trend <- reconstruct(ssa.result, groups = list(Trend = c(1,13:15,18), Seas=c(2:12,16,17,19,20)))
Трендом считаем компоненты- 1, 13, 14, 15, 18; периодикой- 2-12, 16, 17, 19, 20; остальное- шум.
#plot(ssa.rec.trend$Seas)
plot.ts(ts1, type = 'l')
lines(ssa.rec.trend$Trend, col='red')
spectrum(residuals(ssa.rec.trend), log='no', method='pgram', detrend=TRUE)
acf(residuals(ssa.rec.trend))
Тренд выделен хорошо, остаток похож на красный шум, то есть в шуме не осталось тренда и/или периодик.
Модельный ряд, ранг 5 (2 корня от косинуса и 3 от квадратичного тренда).
a <- -1/600
b <- 12
c <- 5
n <- c(1:1000)
ts_mod <- ts(c*exp(a*n)*cos(2/b*pi*n) + (n/200)^2 + 100)
plot(ts_mod, type = 'l')
rank <- 5
ssa <- ssa(ts_mod, L = rank+1, svd.method = 'svd')
r0 <- reconstruct(ssa, groups = list(signal = 1:rank))$signal
plot(ssa)
n_eig <- list(1:rank)
l <- lrr(ssa, group = n_eig)
r1 <- Rssa::roots(l)
r1
## [1] 1.0000252+0.0000429i 1.0000252-0.0000429i 0.9999498+0.0000000i
## [4] 0.8645832+0.4991674i 0.8645832-0.4991674i
plot(l)
parestimate(ssa, group = n_eig)
## period rate | Mod Arg | Re Im
## 145168.017 0.000025 | 1.00003 0.00 | 1.00003 0.00004
## -145168.017 0.000025 | 1.00003 -0.00 | 1.00003 -0.00004
## Inf -0.000051 | 0.99995 0.00 | 0.99995 0.00000
## 12.000 -0.001667 | 0.99833 0.52 | 0.86458 0.49917
## -12.000 -0.001667 | 0.99833 -0.52 | 0.86458 -0.49917
par <- parestimate(ssa, groups = n_eig, method = "esprit")
Полученные 5 корней изобразили на единичной окружности (1 кратности 3, отвечает за квадратичный тренд).
Решая систему получим коэфициенты ЛРФ и нарисуем ряд на фоне модельного.
p <- (par$moduli[1]+par$moduli[2]+par$moduli[3])/3
a <- p
ro <- par$moduli[4]
w <- 2*pi/(par$periods[4])
param <- matrix(c(0, 0, 1, 1, 0,
a, a, a, ro*cos(w), ro*sin(w),
2^2*a^2, 2*a^2, a^2, ro^2*cos(2*w), ro^2*sin(2*w),
3^2*a^3, 3*a^3, a^3, ro^3*cos(3*w), ro^3*sin(3*w),
4^2*a^4, 4*a^4, a^4, ro^4*cos(4*w), ro^4*sin(4*w))
, 5, 5)
c <- solve(t(param), ts_mod[1:5])
c
## [1] 2.499999e-05 4.370836e-05 1.000000e+02 4.322916e+00 -2.495837e+00
t <- ts(Re((c[3]+c[2]*n+c[1]*n^2)*a^n+c[4]*ro^n*cos(w*n)+c[5]*ro^n*sin(w*n)))
plot(t)
lines(ts_mod, col = 'green')
Полученный ряд хорошо описывает исходный.
Модельный ряд с шумом, ранг 5 (2 корня от косинуса и 3 от квадратичного тренда).
ts_mod <- ts_mod + rnorm(1000, 0, 0.15)
plot(ts_mod)
ssa_signal <- ssa(ts_mod, svd.method = 'svd')
r0 <- reconstruct(ssa_signal, groups = list(signal = 1:rank))$signal
plot(ssa_signal, type = 'vectors')
В соответсвии с результатми ssa, ранг = 4, видимо задуманный квадратичный тренд был описан комбинацией двух экспонент.
rank <- 4
n_eig <- list(1:rank)
l1<- lrr(ssa_signal, group = n_eig)
r1 <- Rssa::roots(l1)
main.roots <- r1[Mod(r1) > 0.999]
l2 <- lrr(ssa_signal, reverse = TRUE, group = n_eig)
r2 <- Rssa::roots(l2)
main.roots <- c(main.roots, 1/r2[Mod(r2) > 0.999])
par <- parestimate(ssa_signal, groups = n_eig, method = "esprit")
Отобрали 4 сигнальных корня.
par
## period rate | Mod Arg | Re Im
## Inf 0.000623 | 1.00062 0.00 | 1.00062 0.00000
## Inf -0.000821 | 0.99918 0.00 | 0.99918 0.00000
## 12.000 -0.001682 | 0.99832 0.52 | 0.86457 0.49916
## -12.000 -0.001682 | 0.99832 -0.52 | 0.86457 -0.49916
o <- order(abs(par$periods), decreasing = TRUE)
periods <- (par$periods[o])
moduli <- par$moduli[o]
len <- length(ts_mod)
vars <- matrix(nrow = len, ncol = rank)
for (i in 1:rank) {
if (abs(periods[i]) > 10000)
vars[, i] <- moduli[i]^(1:len)
else if (periods[i] == 2)
vars[, i] <- (-moduli[i])^(1:len)
else if (periods[i] > 0)
vars[, i] <-
moduli[i]^(1:len) * sin(2 * pi * (1:len) / periods[i])
else
vars[, i] <-
moduli[i]^(1:len) * cos(2 * pi * (1:len) / periods[i])
}
lm0 <- lm(r0[1:len] ~ 0 + ., data = data.frame(vars))
coefs0 <- coef(lm0)
coefs0
## X1 X2 X3 X4
## 56.876944333 43.105170602 0.005296981 5.005089155
Нашли коэффициенты ЛРФ, как коэффициенты линейной регрессии.
Нарисуем полученный ряд на фоне модельного.
a1 <- par$moduli[1]
a2 <- par$moduli[2]
ro <- par$moduli[3]
w <- 2*pi/(par$periods[3])
t <- ts(coefs0[1]*a1^n+coefs0[2]*a2^n+coefs0[4]*ro^n*cos(w*n)+coefs0[3]*ro^n*sin(w*n))
plot(t)
lines(ts_mod, col = 'green')
Полученный ряд хорошо описывает исходный.
Применим этот подход к реальному ряду.
plot.ts(ts1, type = 'l')
spectrum(ts1, log='no', method='pgram')
ssa.result <- ssa(ts1, L=336/2)
plot(ssa.result, type = 'vectors', idx = 1:20)
plot(ssa.result, type = 'paired', idx = 2:20)
plot(ssa.result, type = 'series', groups = 1:20)
plot(ssa.result, type = 'series', groups = list(1:20))
plot(ssa.result, type = 'wcor', groups = 1:20)
Ранг равен 15, тренд описывается комбинацией линейной функции и гармоники с большим периодом, сезонность описывается 11 компонентами (пары периодов + “пила”).
rank <- 15
r0 <- reconstruct(ssa.result, groups = list(signal = 1:rank))$signal
n_eig <- list(1:rank)
l1<- lrr(ssa.result, group = n_eig)
r1 <- Rssa::roots(l1)
main.roots <- r1[Mod(r1) > 0.999]
l2 <- lrr(ssa.result, reverse = TRUE, group = n_eig)
r2 <- Rssa::roots(l2)
main.roots <- c(main.roots, 1/r2[Mod(r2) > 0.999])
main.roots
## [1] 1.0029434+0.0000000i -0.5021671+0.8679024i -0.5021671-0.8679024i
## [4] 0.5012392+0.8684362i 0.5012392-0.8684362i 0.8676007+0.5026391i
## [7] 0.8676007-0.5026391i -0.0006864+1.0026209i -0.0006864-1.0026209i
## [10] -0.8678587+0.5016418i -0.8678587-0.5016418i -1.0016702+0.0000000i
plot(main.roots)
par <- parestimate(ssa.result, groups = n_eig, method = "esprit")
par
## period rate | Mod Arg | Re Im
## 6333.541 0.004064 | 1.00407 0.00 | 1.00407 0.00100
## -6333.541 0.004064 | 1.00407 -0.00 | 1.00407 -0.00100
## 2.999 0.002845 | 1.00285 2.10 | -0.50232 0.86797
## -2.999 0.002845 | 1.00285 -2.10 | -0.50232 -0.86797
## 5.999 0.002752 | 1.00276 1.05 | 0.50124 0.86849
## -5.999 0.002752 | 1.00276 -1.05 | 0.50124 -0.86849
## 11.965 0.002649 | 1.00265 0.53 | 0.86756 0.50264
## -11.965 0.002649 | 1.00265 -0.53 | 0.86756 -0.50264
## 2.400 0.002638 | 1.00264 2.62 | -0.86815 0.50160
## -2.400 0.002638 | 1.00264 -2.62 | -0.86815 -0.50160
## 3.999 0.002626 | 1.00263 1.57 | -0.00058 1.00263
## -3.999 0.002626 | 1.00263 -1.57 | -0.00058 -1.00263
## 2.000 0.002112 | 1.00211 3.14 | -1.00211 0.00000
## 92.142 -0.002090 | 0.99791 0.07 | 0.99559 0.06800
## -92.142 -0.002090 | 0.99791 -0.07 | 0.99559 -0.06800
Отобрали сигнальные корни.
p <- (par$moduli[1]+par$moduli[2])/2
o <- order(abs(par$periods), decreasing = TRUE)
periods <- (par$periods[o])
moduli <- par$moduli[o]
len <- length(ts1)
vars <- matrix(nrow = len, ncol = rank)
n <- (1:len)
for (i in 1:rank) {
if (abs(periods[i]) > 1000)
vars[, i] <- n^(i-1)*p^n
else if (periods[i] == 2)
vars[, i] <- (-moduli[i])^n
else if (periods[i] > 0)
vars[, i] <-
moduli[i]^n * sin(2 * pi * n / periods[i])
else
vars[, i] <-
moduli[i]^n * cos(2 * pi * n / periods[i])
}
lm0 <- lm(r0[1:len] ~ 0 + ., data = data.frame(vars))
coefs0 <- coef(lm0)
coefs0
## X1 X2 X3 X4 X5 X6
## 1663.215990 -1.273498 59.028995 41.087106 69.307551 -84.709186
## X7 X8 X9 X10 X11 X12
## 171.204057 -48.993611 155.382530 -39.076783 3.108672 126.328773
## X13 X14 X15
## 9.718763 150.410704 57.620405
Нашли коэффициенты ЛРФ, как коэффициенты линейной регрессии.
Нарисуем полученный ряд на фоне восстановленного.
ro1 <- par$moduli[7]
w1 <- 2*pi/(par$periods[7])
ro2 <- par$moduli[5]
w2 <- 2*pi/(par$periods[5])
ro3 <- par$moduli[11]
w3 <- 2*pi/(par$periods[11])
ro4 <- par$moduli[3]
w4 <- 2*pi/(par$periods[3])
ro5 <- par$moduli[9]
w5 <- 2*pi/(par$periods[9])
ro6 <- par$moduli[13]
w6 <- 2*pi/(par$periods[13])
ro7 <- par$moduli[15]
w7 <- 2*pi/(par$periods[15])
t <- ts((coefs0[1]+coefs0[2]*n)*p^n+
coefs0[3]*ro7^n*cos(w7*n)+coefs0[4]*ro7^n*sin(w7*n)+
coefs0[5]*ro1^n*cos(w1*n)+coefs0[6]*ro1^n*sin(w1*n)+
coefs0[7]*ro2^n*cos(w2*n)+coefs0[8]*ro2^n*sin(w2*n)+
coefs0[9]*ro3^n*cos(w3*n)+coefs0[10]*ro3^n*sin(w3*n)+
coefs0[11]*ro4^n*sin(w4*n)+coefs0[12]*ro4^n*cos(w4*n)+
coefs0[13]*ro5^n*sin(w5*n)+coefs0[14]*ro5^n*cos(w5*n)+
coefs0[15]*ro6^n*cos(w6*n)
, start = c(1992, 1), frequency = 12)
plot(r0)
lines(t, col = 'green')
plot(ts1)
lines(t, col = 'red')
Полученный ряд хорошо описывает восстановленный ряд.
Для начала зададим число групп, первый вариант– 4 групп (1 отвечает за тренд и 3 за сезонность), второй вариант– 6 групп (1 отвечает за тренд и 5 за сезонность),
plot(ssa.result, type = "vectors", idx = 1:20)
plot(ssa.result, type = "series", groups = 1:20)
plot(ssa.result, type = "wcor", groups = 1:20)
g <- grouping.auto(ssa.result, freq.bins = 4, threshold = 0)
plot(reconstruct(ssa.result, groups = g))
g <- grouping.auto(ssa.result, freq.bins = 6, threshold = 0)
plot(reconstruct(ssa.result, groups = g))
В первом случае в тренд попала сезонность, значит количество групп надо увеличить, во втором случае тренд выделен хорошо.
Второй случай– зададим freq.bins = 1/(12*2) (меньше 1/12- самой маленькой возможной частоты из сезонности) и threshold = 0 для определения порога.
g <- grouping.auto(ssa.result, freq.bins = list(Trend = 1/(12*2)), threshold = 0)
plot(reconstruct(ssa.result, groups = g))
plot(g, type='b', order = TRUE)
Будем считать компоненты частью тренда при привышении порога в 0.85.
g <- grouping.auto(ssa.result, freq.bins = list(Trend = 1/(12*2)), threshold = 0.85)
plot(ts1)
lines(reconstruct(ssa.result, groups = g)$Trend, col='red')
g$Trend
## [1] 1 13 14 15 18 31 32
К тренду отнесены компоненты 1, 13, 14, 15, 18, 31, 32. Тренд выделен хорошо.
fossa помогает, если есть слабая разделимость, но нет сильной. Здесь смешаны компоненты 13-15.
f <- fossa(ssa.result, nested.group = 13:15)
plot(ssa.result, type = "vectors", idx = 1:15)
plot(f, type = "vectors", idx = 1:15)
plot(ssa.result, type = "series", groups = 1:15)
plot(f, type = "series", groups = 1:15)
plot(ssa.result, type = "wcor", groups = 1:15)
plot(f, type = "wcor", groups = 1:15)
Не помогло :(
iossa уже и с отсутствием слабой разделимости может справиться. Подаем на вход компоненты, отнесенные к тренду и компоненты, отнесенные к периодике.
g <- list(c(1, 13:15), 2:12)
issa <- iossa(ssa.result, nested.group = g)
plot(ssa.result, type = "vectors", idx = 1:15)
plot(issa, type = "vectors", idx = 1:15)
## Warning in .contribution(x, idx, ...): Elementary matrices are not F-orthogonal
## (max F-cor is -0.000776). Contributions can be irrelevant
plot(ssa.result, type = "series", groups = 1:15)
plot(issa, type = "series", groups = 1:15)
r0 <- reconstruct(issa, groups = issa$iossa.groups)
plot(r0)
Собственные вектора 13-15 сгладились.
eossa также справляется с отсутствием сильной разделимости, на вход подаем количество групп = 8.
e <- eossa(ssa.result, nested.group = g, k = 8)
plot(ssa.result, type = "vectors", idx = 1:15)
plot(e, type = "vectors", idx = 1:15)
## Warning in .contribution(x, idx, ...): Elementary matrices are not F-orthogonal
## (max F-cor is -0.0119). Contributions can be irrelevant
plot(e, type = "series", groups = e$iossa.groups)
r0 <- reconstruct(e, groups = e$iossa.groups)
plot(r0)
Собственные вектора 13-15 также сгладились, в группы были объедены компоненты, отвечающие за линейную составлющую тренда, за гармонику с большим периодом из тренда и все периодики из сезонности попарно.
Рассмотрим два варианта: “отрежем” один и два периода и будем прогнозировать по оставшимся точкам “отрезанные”.
ts1s <- ts(ts1[1:324], start = c(1992, 1), frequency = 12)
ts1s2 <- ts(ts1[1:312], start = c(1992, 1), frequency = 12)
ssa.result_1 <- ssa(ts1s, L=324/2)
ssa.result_2 <- ssa(ts1s2, L=312/2)
plot(ssa.result_1, type = 'vectors', idx = 1:20)
plot(ssa.result_1, type = 'paired', idx = 2:20)
plot(ssa.result_1, type = 'series', groups = 1:20)
plot(ssa.result_1, type = 'series', groups = list(1:20))
plot(ssa.result_1, type = 'wcor', groups = 1:20)
plot(ssa.result_2, type = 'vectors', idx = 1:20)
plot(ssa.result_2, type = 'paired', idx = 2:20)
plot(ssa.result_2, type = 'series', groups = 1:20)
plot(ssa.result_2, type = 'series', groups = list(1:20))
plot(ssa.result_2, type = 'wcor', groups = 1:20)
При “отрезании” одного периода при группировке будем рассматривать 15 компонент.
При “отрезании” двух периодов при группировке будем рассматривать 12 компонент, так как дальше компоненты сильно смешаны.
r <- rforecast(ssa.result_1, groups = list(1:15), len = 12)
plot(ts(c(ts1s,r), start=start(ts1), frequency=frequency(ts1)))
lines(ts1, col = 'red')
rec_1 <- rmse(ts1[325:336], r)
rec_1
## [1] 108.9153
r <- rforecast(ssa.result_2, groups = list(1:12), len = 24)
plot(ts(c(ts1s2,r), start=start(ts1), frequency=frequency(ts1)))
lines(ts1, col = 'red')
rec_2 <- rmse(ts1[313:336], r)
rec_2
## [1] 134.0593
Ошибка рекуррентного прогноза для 1 периода– 108.92, для 2 периодов– 134.06.
v <- vforecast(ssa.result_1, groups = list(1:15), len = 12)
plot(ts(c(ts1s,v), start=start(ts1), frequency=frequency(ts1)))
lines(ts1, col = 'red')
vec_1 <- rmse(ts1[325:336], v)
vec_1
## [1] 85.11311
v <- vforecast(ssa.result_2, groups = list(1:12), len = 24)
plot(ts(c(ts1s2,v), start=start(ts1), frequency=frequency(ts1)))
lines(ts1, col = 'red')
vec_2 <- rmse(ts1[313:336], v)
vec_2
## [1] 134.5759
Ошибка векторного прогноза для 1 периода– 85.11, для 2 периодов– 134.58.
bf_r <- forecast(ssa.result_1, groups = list(1:15), method = "recurrent", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_r)
lines(ts1, col = 'black', t='l')
brec_1 <- rmse(ts1[325:336], bf_r$mean)
brec_1
## [1] 108.9153
bf_v <- forecast(ssa.result_1, groups = list(1:15), method = "vector", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_v)
lines(ts1, col = 'black', t='l')
bvec_1 <- rmse(ts1[325:336], bf_v$mean)
bvec_1
## [1] 85.11311
Ошибка bootstrap прогноза для 1 периода:
Рекуррентный– 108.92.
Векторный– 85.11.
Обычные и рекуррентный и векторный прогнозы показали такие же результаты.
Также построены 95%-предсказательные интервалы для обоих прогнозов.
bf_r <- forecast(ssa.result_2, groups = list(1:12), method = "recurrent", bootstrap = TRUE, len = 24, R = 100, level = 0.95, interval = 'prediction')
plot(bf_r)
lines(ts1, col = 'black', t='l')
brec_2 <- rmse(ts1[313:336], bf_r$mean)
brec_2
## [1] 134.0593
bf_v <- forecast(ssa.result_2, groups = list(1:12), method = "vector", bootstrap = TRUE, len = 24, R = 100, level = 0.95, interval = 'prediction')
plot(bf_v)
lines(ts1, col = 'black', t='l')
bvec_2 <- rmse(ts1[313:336], bf_v$mean)
bvec_2
## [1] 134.5759
Ошибка bootstrap прогноза для 2 периода:
Рекуррентный– 134.06.
Векторный– 134.58.
Нет особых различий между векторным и рекуррентным прогнозами.
Обычные и рекуррентный и векторный прогнозы показали такие же результаты.
Также построены 95%-предсказательные интервалы для обоих прогнозов.
На примере для двух периода:
Плохо отделяется тренд, для исправления этого используем iossa.
g <- list(c(1, 13:15), 2:12)
issa <- iossa(ssa.result_2, nested.group = g)
plot(ssa.result_2, type = "vectors", idx = 1:15)
plot(issa, type = "vectors", idx = 1:15)
## Warning in .contribution(x, idx, ...): Elementary matrices are not F-orthogonal
## (max F-cor is -0.00102). Contributions can be irrelevant
plot(ssa.result_2, type = "series", groups = 1:15)
plot(issa, type = "series", groups = 1:15)
Первые 4 собственных вектора, соответствующих компонентам тренда стали более гладкие.
Сравним прогноз тренда до и после.
Для сравнения выделим тренд каким-то другим способом на полных данных (рассмотрим ранее полученные результаты при помощи ssa на полных данных).
До:
bf_r <- forecast(ssa.result_2, groups = list(c(1, 13:15)), method = "recurrent", bootstrap = TRUE, len = 24, R = 100, level = 0.95)
plot(bf_r)
lines(ts1, col = 'black', t='l')
lines(ssa.rec.trend$Trend, col='red', t='l')
brec_2 <- rmse(ssa.rec.trend$Trend[313:336], bf_r$mean)
brec_2
## [1] 125.9559
bf_v <- forecast(ssa.result_2, groups = list(c(1, 13:15)), method = "vector", bootstrap = TRUE, len = 24, R = 100, level = 0.95)
plot(bf_v)
lines(ts1, col = 'black', t='l')
lines(ssa.rec.trend$Trend, col='red')
bvec_2 <- rmse(ssa.rec.trend$Trend[313:336], bf_v$mean)
bvec_2
## [1] 80.70864
После:
bf_r <- forecast(issa, groups = list(1:4), method = "recurrent", bootstrap = TRUE, len = 24, R = 100, level = 0.95)
plot(bf_r)
lines(ts1, col = 'black', t='l')
lines(ssa.rec.trend$Trend, col='red')
brec_2 <- rmse(ssa.rec.trend$Trend[313:336], bf_r$mean)
brec_2
## [1] 172.6782
bf_v <- forecast(issa, groups = list(1:4), method = "vector", bootstrap = TRUE, len = 24, R = 100, level = 0.95)
plot(bf_v)
lines(ssa.rec.trend$Trend, col='red')
bvec_2 <- rmse(ssa.rec.trend$Trend[313:336], bf_v$mean)
bvec_2
## [1] 133.5523
Стало хуже :)
На восстановленных рядах 5, 7, 11, 13 виден резкий скачок в конце, а полученный тренд как раз находится значительно ниже тренда, полученного ssa на полных данных. Возможно причина именно в этом.
rank <- 15
r0 <- reconstruct(ssa.result_1, groups = list(signal = 1:rank))$signal
n_eig <- list(1:rank)
l1<- lrr(ssa.result_1, group = n_eig)
r1 <- Rssa::roots(l1)
main.roots <- r1[Mod(r1) > 0.999]
l2 <- lrr(ssa.result_1, reverse = TRUE, group = n_eig)
r2 <- Rssa::roots(l2)
main.roots <- c(main.roots, 1/r2[Mod(r2) > 0.999])
main.roots
## [1] -0.5021028+0.8682899i -0.5021028-0.8682899i 1.0028899+0.0000000i
## [4] 0.8676498+0.5026699i 0.8676498-0.5026699i 0.5011646+0.8684447i
## [7] 0.5011646-0.8684447i -0.8681241+0.5014793i -0.8681241-0.5014793i
## [10] -0.0006135+1.0023935i -0.0006135-1.0023935i -1.0016530+0.0000000i
par <- parestimate(ssa.result_1, groups = n_eig, method = "esprit")
par
## period rate | Mod Arg | Re Im
## Inf 0.006012 | 1.00603 -0.00 | 1.00603 -0.00000
## Inf 0.003966 | 1.00397 0.00 | 1.00397 0.00000
## 2.999 0.003149 | 1.00315 2.10 | -0.50229 0.86834
## -2.999 0.003149 | 1.00315 -2.10 | -0.50229 -0.86834
## 2.400 0.002799 | 1.00280 2.62 | -0.86841 0.50147
## -2.400 0.002799 | 1.00280 -2.62 | -0.86841 -0.50147
## 11.969 0.002789 | 1.00279 0.52 | 0.86776 0.50258
## -11.969 0.002789 | 1.00279 -0.52 | 0.86776 -0.50258
## 5.999 0.002781 | 1.00278 1.05 | 0.50120 0.86855
## -5.999 0.002781 | 1.00278 -1.05 | 0.50120 -0.86855
## 3.999 0.002534 | 1.00254 1.57 | -0.00048 1.00254
## -3.999 0.002534 | 1.00254 -1.57 | -0.00048 -1.00254
## 2.000 0.002153 | 1.00216 3.14 | -1.00216 0.00000
## 92.145 -0.001857 | 0.99814 0.07 | 0.99582 0.06801
## -92.145 -0.001857 | 0.99814 -0.07 | 0.99582 -0.06801
Тренд описывается четырьмя компонентами; корень, соответствующий линейному тренду имеет модуль больше 1, значит линейный тренд возрастает, а корень, соответствующий гармонике, относящейся к тренду имеет модуль меньше 1, значит колебания затухают. Все корни, соответствующие гармоникам имеют модули больше 1, а значит колебания усиливаются.
Проверим это, спрогнозировав на 20 периодов вперед:
bf_v <- forecast(ssa.result_1, groups = list(1:15), method = "vector", bootstrap = TRUE, len = 12*20, R = 100)
plot(bf_v)
Тренд возрастает, а колебания увеличиваются.
rank <- 12
r0 <- reconstruct(ssa.result_2, groups = list(signal = 1:rank))$signal
n_eig <- list(1:rank)
l1<- lrr(ssa.result_2, group = n_eig)
r1 <- Rssa::roots(l1)
main.roots <- r1[Mod(r1) > 0.999]
l2 <- lrr(ssa.result_2, reverse = TRUE, group = n_eig)
r2 <- Rssa::roots(l2)
main.roots <- c(main.roots, 1/r2[Mod(r2) > 0.999])
main.roots
## [1] 1.0031701+0.0000000i -0.5017732+0.8682682i -0.5017732-0.8682682i
## [4] -0.8684850+0.5012696i -0.8684850-0.5012696i 0.8675172+0.5026221i
## [7] 0.8675172-0.5026221i 0.5012153+0.8682713i 0.5012153-0.8682713i
## [10] -0.0004278+1.0024724i -0.0004278-1.0024724i -1.0022579+0.0000000i
par <- parestimate(ssa.result_2, groups = n_eig, method = "esprit")
par
## period rate | Mod Arg | Re Im
## Inf 0.003225 | 1.00323 0.00 | 1.00323 0.00000
## 2.999 0.003017 | 1.00302 2.09 | -0.50199 0.86837
## -2.999 0.003017 | 1.00302 -2.09 | -0.50199 -0.86837
## 2.400 0.002904 | 1.00291 2.62 | -0.86864 0.50129
## -2.400 0.002904 | 1.00291 -2.62 | -0.86864 -0.50129
## 11.971 0.002738 | 1.00274 0.52 | 0.86777 0.50246
## -11.971 0.002738 | 1.00274 -0.52 | 0.86777 -0.50246
## 5.999 0.002707 | 1.00271 1.05 | 0.50125 0.86843
## -5.999 0.002707 | 1.00271 -1.05 | 0.50125 -0.86843
## 2.000 0.002630 | 1.00263 3.14 | -1.00263 0.00000
## 3.999 0.002615 | 1.00262 1.57 | -0.00037 1.00262
## -3.999 0.002615 | 1.00262 -1.57 | -0.00037 -1.00262
Тренд описывается одной компонентой, соответствующий корень имеет модуль больше 1, значит тренд возрастающий. Все корни, соответствующие гармоникам также имеют модули больше 1, а значит колебания усиливаются.
Проверим это, спрогнозировав на 20 периодов вперед:
bf_v <- forecast(ssa.result_2, groups = list(1:12), method = "vector", bootstrap = TRUE, len = 12*20, R = 100)
plot(bf_v)
Тренд возрастает, а колебания увеличиваются.
Генерируем стационарный ряд + шум.
b <- 12
n <- 1:120
ts_mod <- ts(cos(2/b*pi*n) + 2 + rnorm(120,0,0.07))
plot(ts_mod)
ssa_toep1 <- ssa(ts_mod, L = 60, kind = '1d-ssa')
ssa_toep2 <- ssa(ts_mod, L = 60, kind = 'toeplitz-ssa')
plot(ssa_toep1, type = 'vectors', idx = 1:5)
plot(ssa_toep1, type = 'paired', idx = 2:5)
plot(ssa_toep1, type = 'series', groups = 1:5)
plot(ssa_toep1, type = 'series', groups = list(5))
plot(ssa_toep1, type = 'wcor', groups = 1:5)
plot(ssa_toep2, type = 'vectors', idx = 1:5)
plot(ssa_toep2, type = 'paired', idx = 2:5)
plot(ssa_toep2, type = 'series', groups = 1:5)
plot(ssa_toep2, type = 'series', groups = list(1:5))
plot(ssa_toep2, type = 'wcor', groups = 1:5)
Сигнальными компонентами считаем первые 3– первая компонента- тренд, вторая и третья- косинус.
ssa.rec_toep1 <- reconstruct(ssa_toep1, groups = list(Trend = 1, Seas=c(2:3)))
ssa.rec_toep2 <- reconstruct(ssa_toep2, groups = list(Trend = 1, Seas=c(2:3)))
Посмотрим насколько восстановленный ряд отличается от модельного ряда.
ssa_signal1 <- ssa.rec_toep1$Trend + ssa.rec_toep1$Seas
ssa_signal2 <- ssa.rec_toep2$Trend + ssa.rec_toep2$Seas
rmse(ts_mod, ssa_signal1)
## [1] 0.06722906
rmse(ts_mod, ssa_signal2)
## [1] 0.09872395
Basic-SSA показывает результаты лучше, чем Toeplitz-SSA, уже при таком небольшом уровне шума.
ssa.res_2center <- ssa(ts1, L=336/2, row.projector = "center", column.projector = "center")
plot(ssa.res_2center, type = 'vectors', idx = 1:25)
plot(ssa.res_2center, type = 'paired', idx = 3:25)
plot(ssa.res_2center, type = 'series', groups = 1:25)
plot(ssa.res_2center, type = 'series', groups = list(1:25))
plot(ssa.res_2center, type = 'wcor', groups = 1:25)
Трендом считаем компоненты– 1, 2, 14, 15, 16, 19, 22; сезонностью– 3-13, 17, 18, 20, 21.
ssa.rec <- reconstruct(ssa.res_2center, groups = list(Trend = c(1,2,14:16,19,22), Seas=c(3:13,17:18,20:21)))
plot.ts(ts1, type = 'l')
lines(ssa.rec$Trend, col='red')
lines(ssa.rec.trend$Trend, col='green')
spectrum(residuals(ssa.rec), log='no', method='pgram', detrend=TRUE)
acf(residuals(ssa.rec))
Сравним тренд, выделеный классическим ssa и ssa с двойным центрированием: на концах заметны отличия, адекватнее тренд описывает классический ssa, возможно, это связано с тем, что данный ряд имеет нелинейный тренд, так как в нем присутствует гармоника с большим периодом.
Остаток похож на красный шум, но для этого пришлось брать большее число компонент, чем при классическом ssa.
Добавим пропуски.
ts_gap <- ts1
gap <- 193:246
ts_gap[gap] <- NA
plot(ts_gap)
Рассмотрим длину окна такую что, вторая часть имеет не ноль векторов вложений.
ssa_gap_res <- ssa(ts_gap, L=72)
plot(ssa_gap_res, type = 'vectors', idx = 1:20)
plot(ssa_gap_res, type = 'paired', idx = 2:20)
plot(ssa_gap_res, type = 'series', groups = 1:20)
plot(ssa_gap_res, type = 'series', groups = list(1:20))
plot(ssa_gap_res, type = 'wcor', groups = 1:20)
Результат применения gapfill:
gr <- list(c(1:12, 13, 14, 19, 20))
g <- gapfill(ssa_gap_res, groups = gr, method = "sequential",
base = "reconstructed")
plot(ts_gap, col = "black")
lines(g, col = "red")
plot(ts1, col = "black")
lines(ts(g[gap], start = c(2008, 1), frequency = 12), col = "red")
rmse(g[gap], ts1[gap])
## [1] 99.42902
ig <- igapfill(ssa_gap_res, groups = gr,
base = "reconstructed", maxiter = 2000)
plot(ts_gap, col = "black")
lines(ig, col = "red")
plot(ts1, col = "black")
lines(ts(ig[gap], start = c(2008, 1), frequency = 12), col = "red")
rmse(ig[gap], ts1[gap])
## [1] 93.25305
igapfill c maxiter=2000 показывает результат лучше (с точки зрения rmse), чем gapfill.
Аппроксимируем исходный ряд рядом конечного ранга:
cadz <- cadzow(ssa.result, rank = 15, tol = 1e-10)
plot(ts1)
lines(cadz, col = 'red')
rmse(ts1, cadz)
## [1] 79.63702
plot(ssa(cadz))
Получили ряд конечного ранга = 15.
Рассмотрим вариант с весами, для достижения большей точности.
L <- 336/2
K <- length(ts1) - L + 1
alpha <- 0.01
weights <- vector(len = K)
weights[1:K] <- alpha
weights[seq(K, 1, -L)] <- 1
ssa_cadz_weights <- ssa(ts1, L = L, column.oblique = "identity",
row.oblique = weights)
cadz_weights <- cadzow(ssa_cadz_weights, rank = 15, maxiter = 10)
plot(ts1)
lines(cadz_weights, col = 'red')
rmse(ts1, cadz_weights)
## [1] 68.53668
plot(ssa(cadz_weights))
rmse меньше по сравнению с обычным Cadzow.
K.sliding <- 2
forecast.base.len <- 2*frequency(ts1s)
base.len <- length(ts1s)
sliding.len <- base.len - K.sliding - forecast.base.len + 1
ncomp <- 1:22
L.min <- 4*frequency(ts1s)
Ls <- seq(L.min, length(ts1s)/2, by = frequency(ts1s))
m0 <- forecast.sliding.rmse(ts1s,
K.sliding = K.sliding,
L = Ls, ncomp = ncomp,
method = "recurrent",
forecast.len = forecast.base.len,
.progress = "none")
par <- optim.par(m0)
par$L.opt
## [1] 96
par$ncomp.opt
## [1] 19
par$m
##
## L 1 2 3 4 5 6 7 8
## 48 299953.1 233286.9 209547.6 171339.2 147837.7 109897.4 104276.8 84959.12
## 60 300770.9 234748.7 209677.6 170746.5 147855.6 109503.2 104222.1 82925.58
## 72 302237.5 236795.8 210847.1 171382.3 149042.5 110470.1 105539.2 82910.69
## 84 304244.0 239164.3 212749.6 173051.5 150872.3 112173.2 107477.1 84272.42
## 96 306387.1 241653.4 214909.2 175252.9 153017.0 114315.6 109536.0 86265.41
## 108 307353.6 242671.5 215892.9 176053.8 153983.8 115229.7 110586.8 87142.54
## 120 307666.1 243108.7 216192.7 176229.1 154257.1 115393.7 111138.9 87432.56
## 132 308207.8 243738.7 216734.6 176719.8 154825.6 115917.6 111898.0 87566.09
## 144 309398.5 244830.0 218041.6 178012.8 156068.5 117283.2 113231.9 88410.93
## 156 309888.9 245423.7 218581.7 178471.3 156625.0 117575.2 113664.7 88719.02
##
## L 9 10 11 12 13 14 15
## 48 49754.03 37248.92 18130.66 10231.960 9687.210 6494.669 6028.608
## 60 48578.32 33844.42 16511.54 9067.978 8032.647 5893.959 7281.058
## 72 49511.84 32891.48 17099.49 9871.432 7895.322 8241.512 8317.616
## 84 51214.10 33827.41 18672.66 11617.241 7961.980 9283.513 7855.896
## 96 53270.65 35385.83 20812.33 13823.281 8583.109 13106.236 12130.460
## 108 54502.91 36572.21 22129.72 15225.695 9593.236 17043.091 14907.630
## 120 55160.41 37004.88 22963.98 16007.474 12661.867 19583.756 26913.524
## 132 55860.01 37605.57 23905.51 16866.262 15618.963 25711.432 38682.329
## 144 56987.76 39293.25 25145.77 18147.784 19164.729 27267.362 31318.891
## 156 57173.76 39604.56 25442.52 18487.078 22806.867 26377.526 32160.950
##
## L 16 17 18 19 20 21 22
## 48 9539.674 8417.259 10341.970 10267.722 10008.970 11075.830 12608.131
## 60 7585.594 6416.566 6662.913 7452.579 8225.342 7565.183 12687.272
## 72 6158.725 5421.843 5165.162 6876.015 6794.874 6636.124 6226.282
## 84 7319.978 6057.928 4984.602 4840.112 6632.341 5832.768 5442.083
## 96 11265.733 5398.528 4362.144 4176.981 11777.094 10947.351 10365.682
## 108 14219.807 9028.584 8169.678 6615.222 4930.718 5243.344 5939.900
## 120 21554.983 19104.077 17760.675 16135.923 4229.633 4930.668 8401.067
## 132 36407.412 29051.667 27696.440 27311.061 15029.839 12463.532 40713.008
## 144 28135.336 24524.378 22608.222 21664.438 4676.953 94202.460 173117.126
## 156 29603.955 24596.428 21956.432 21456.467 7286.060 67050.744 146725.683
matplot(Ls, sqrt(par$m), ylab = "", xlab = "Window lengths",
type = "l", col = topo.colors(10))
Построим прогноз (векторный и рекуррентный bootstrap) с оптимальными параметрами:
ssa.result_optim <- ssa(ts1s, L=96)
plot(ssa.result_optim, type = 'vectors', idx = 1:19)
plot(ssa.result_optim, type = 'paired', idx = 2:19)
plot(ssa.result_optim, type = 'series', groups = 1:19)
plot(ssa.result_optim, type = 'series', groups = list(1:19))
plot(ssa.result_optim, type = 'wcor', groups = 1:19)
bf_r <- forecast(ssa.result_optim, groups = list(1:19), method = "recurrent", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_r)
lines(ts1, col = 'black', t='l')
brec_3 <- rmse(ts1[325:336], bf_r$mean)
brec_3
## [1] 74.5293
bf_v <- forecast(ssa.result_optim, groups = list(1:19), method = "vector", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_v)
lines(ts1, col = 'black', t='l')
bvec_3 <- rmse(ts1[325:336], bf_v$mean)
bvec_3
## [1] 72.34684
rmse векторного прогноза немного меньше.
Разладка есть в данных (последние два периода).
plot(ts)
N <- length(ts)
ssa_full_res <- ssa(ts, L = N/2)
plot(ssa_full_res, type = 'vectors', idx = 1:20)
plot(ssa_full_res, type = 'paired', idx = 2:20)
plot(ssa_full_res, type = 'series', groups = 1:20)
plot(ssa_full_res, type = 'series', groups = list(1:20))
plot(ssa_full_res, type = 'wcor', groups = 1:20)
Рассмотрим первые 12 компонент:
rank <- 12
M <- 120
L <- M/4
hm <- hmatr(ts, B = M, T = 36, L = L, neig = rank)
plot(hm)
На матрице неоднородности видна разладка в последних периодах.